Emergence of Hopf bifurcation in an extended SIR dynamic

In this paper, the original SIR model is improved by considering a new compartment, representing the hospitalization of critical cases. A system of differential equations with four blocks is developed to analyze the treatment of severe cases in an Intensive Care Unit (ICU). The outgoing rate of the infected individuals who survive is divided into nI and bII+b where the second term represents the transition rate of critical cases that are hospitalized in ICU. The findings demonstrate the existence of forward, backward and Hopf bifurcations in various ranges of parameters.


Introduction
Mckendrick [1] developed the first Susceptible-Infectious-Recovered (SIR) model to simulate and predict a disease spread phenomenon and its epidemic state. In this standard model, there are three blocks labeled S, I and R which represent the number or percentage of susceptible, infectious and recovered individuals. The transitions between these compartments are happened by constants β and γ, respectively indicating the infection rate and the recovery rate.
Since the development of the original SIR model, the parameters β and γ have been modified in various studies. Some new blocks have been introduced in order to develop better models to more accurately predict the behavior of different epidemics. For instance, β was replaced with kI p−1 S q−1 , β aS+bI+cR , μe −mI in [2][3][4], respectively where the researchers developed the model of individuals' behavior along with government strategies by considering β as a function of other blocks. More recent models on COVID-19 can be found in the articles [5][6][7].
Moreover, γ has been altered in different studies to consider the capacity of the general healthcare system in a country. For instance, γI was replaced with: in the articles [8][9][10], respectively. In our study, an improved SIR model that examines the effects of ICU in hospitals was deployed. The importance of Hospitalization and ICU in real situations can be found also in the articles [11][12][13]. To this end, a new block H was introduced to express the number of individuals in ICU. There is, hence, a new block in the differential equations of the SIR model. Fig  1 represents the schematic of this model with its differential equations expressed as follows: In the Eq (4), A is the birth rate, d is the natural death rate, α is the death rate caused by a disease, and β is the incidence rate. In this proposed model, individuals of I who survive are divided into two groups. The individuals of the first group are cured without hospitalization and do not need special medical equipment utilized in an ICU. This group of patients are transported to the block R with the transition rate nI. On the other hand, the subjects in the second group need to receive special medical treatments in hospital. These individuals are  Tables 1 and 2. https://doi.org/10.1371/journal.pone.0276969.g001 transferred to a new block H with a transition rate μ(I). The individuals of H are recovered with the rate n 0 H. remark 0.1. Obviously, the parameters n and n 0 can be functions of I and H but in order to simplify the equations and future analysis, they are considered as constant. Furthermore, α, which is considered as constant like previous argument, represents all the death rate caused by infectious for both of normal and severe patients. One could introduce a different death rate for severe cases which does not actually affect the behavior of the system as we show in future theorems. Therefore, it is more reasonable and convenient to consider just one (average) death rate caused by infectious in block I and to concentrate on the transition function to the block H, μ(I), which models the impact of sever patients on behavior of the system and represents the interaction between normal and sever cases.
A function which can model μ(I) in real situations can be expressed as follows, where b represents the number of beds in ICU: At first, when the number of infected individuals I is less than the number of beds b, the transition rate to ICU is I. For example, there are 100 beds and 50 infected people, so exactly 50 individuals are transported to ICU. However, when the number of infected individuals I exceeds the number of beds b, the transition rate to ICU is independent of the number of infected individuals and is identified as b. For example, there are 150 infected people and 100 beds, so 100 infected people are transferred to ICU.
However, this procedure has a serious shortcoming, as the function is obviously not differentiable at the point I = b. Thus the following function was selected which is both similar to Obviously, in the limits I ! 0 and I ! 1, we have μ 2 (I)!μ 1 (I) and dm 2 dI ! dm 1 dI . Therefore, μ (I) is defined as follows: The differential equations are expressed as follows: In the system (5), the first two equations are independent of the other equations, so just these equations are considered in order to analyze the system (5) and find its fixed points and bifurcations. As such, the essential equations are: where δ is defined as follows: This article consists of just mathematical analysis and all of the theorems are proved by mathematical analysis and not by numerical simulations. Furthermore, the aim of this article is to explain and research local behavior of the introduced model and does not explain its applications to the real situations with help of numerical analysis. However, there are several numerical simulations related to some proved theorems in order to clarify their possible applications and to give intuition, and not to justify them.

Fixed points
If we consider the system (6) as _ X ¼ f ðXÞ, we must solve f(X � ) = 0 in order to find fixed points, then: By solving the second equation, one can find: Next by combining the first case and the first equation, one fixed point can be determined as follows: By combining the second case and the first equation, the following quadratic equation can be formulated: As it will be discussed in the following sections, it is advisable to define R 0 , basic reproduction number, as follows: Now with this definition, the Eq (8) can be rewritten as follows: with the following coefficients: As a result, the following theorems can be formulated about other fixed points X � 6 ¼ E 0 of the system (6).
Proof. It is obvious that: Therefore, the quadratic Eq (10) has two solutions I 1 and I 2 with the condition I 1 I 2 < 0. However, the variable I must be non-negative in the system (6), so there is just one permissible solution for the quadratic Eq (10) and this proves the theorem.
It is obvious that R 0 < 1 so c 0 > 0, and we can infer from the assumption that: On the one hand, if there are real solutions for the quadratic Eq (10), both of them must be negative (c 1 > 0 & c 0 > 0). On the other hand, there are no real solutions for this equation. Therefore, the theorem has been proved.
Theorem 0.3. If d ðdþ1Þ < R 0 < 1, there are the following cases: : there are two fixed points X � when R 0 is close enough to 1.
Proof. If R 0 < 1, it is obvious for the two cases that 0 < c 0 . Now we can derive the following result for the first case: As c 1 > 0 and c 0 > 0, the proof of this case is similar to the previous theorem. Thus there is no fixed point X � in this case.
The second case is somewhat complicated. At first, the new parameter � is defined as follows: � is positive, because: 0 < � and the following properties are obvious: However, as b < d bðdþ1Þ in the second case, the right side of the last inequality is positive, so � can be found. As a result: By considering the properties of �, which is mentioned above, we can find that the right inequality is equal to the condition, with R 0 close enough to 1.
When c 1 < 0 and c 2 > 0, I min ≔ argmin P(I) is positive and when R 0 = 1, there are two solutions for the quadratic Eq (10): I 1 = 0 and I 2 > 0. Therefore, when R 0 = 1, P(I min )<0. Now by decreasing R 0 , when R 0 is close enough to 1 and the conditions c 1 < 0 and P(I min )<0 are held, there must be two positive solutions for the quadratic Eq (10). (The condition P(I min )<0 can be held because of continuity.) In this way, the second case of the theorem can be proved.

Stability of fixed points
First, we must determine Df in order to analyze the stability of fixed points. Df is expressed as follows: Now Df(E 0 ) can be expressed as follows: Therefore, the following theorem is resulted: is an upper triangular matrix, its eigenvalues are −d and (δ + 1)(R 0 − 1). Now if R 0 < 1, then there are two negative eigenvalues and E 0 will be a stable node; and if R 0 > 1, then there is one negative and one positive eigenvalue and E 0 will be a saddle point.
When other fixed points X � exist, Df(X � ) can be expressed as follows: Now the trace of Df(X � ), which is the sum of the eigenvalues of Df(X � ), is expressed as follows: The determinant of Df(X � ), which is the product of the eigenvalues of Df(X � ), can be expressed by the following equations: There are the following theorems about the stability of fixed points X � : Theorem 0.5. When R 0 > 1, the unique fixed point X � is stable if 1 4 Proof. First, the necessary condition, which must be satisfied in both of the real and complex eigenvalues cases, is λ 1 + λ 2 < 0, i.e: It is now easy to realize that the function Ib ðIþbÞ 2 is less than or equal to 1 4 for non-negative I, and as we saw in the theorem, there is a fixed point with I 2 (0, +1) when R 0 > 1. Therefore, the above inequality must be satisfied for 8I 2 (0, +1), so: However, when the eigenvalues of Df(X � ) are real, another condition that must be satisfied is λ 1 λ 2 > 0, i.e: One can observe that the above inequality is satisfied in this interval if and only if: And this completes the proof of this theorem.

Bifurcations
In this section, it is shown that the system (6) undergoes three kinds of bifurcation, forward and backward bifurcations and Hopf bifurcation.

Forward & backward bifurcation
Theorem 0.6 (forward bifurcation). When R 0 = 1, the system (6) undergoes a forward bifurcation if 1 4 < d and d ðdþ1Þb < b. Proof. According to the theorems of the previous sections, when d ðdþ1Þb < b, E 0 is a unique stable fixed point for R 0 < 1 and when 1 4 < d and d ðdþ1Þb < b, there is a unique stable fixed point X � for R 0 > 1 and E 0 is an unstable fixed point. Therefore, when d ðdþ1Þb < b and 1 4 < d, both of the conditions are satisfied and there is a forward bifurcation in the system (6).
As an example, this theorem is demonstrated in the plots of Fig 2. Theorem 0.7 (backward bifurcation). The system (6) undergoes a backward bifurcation for Proof. According to the theorems in the previous sections, when R 0 < 1, E 0 is always a stable fixed point and there are two other fixed points X � if b < d ðdþ1Þb and R 0 close enough to 1. Now we claim that one of these fixed points X � is a stable node and the other one is a saddle point. We prove this by considering the Index Theory. First, we choose one closed curve C like in Fig 3: This curve C is chosen large enough, so when the fixed points X � exist, this curve C encloses them. (The semicircle in the fourth quadrant is too small, so it encloses just E 0 for R 0 < 1). Now when R 0 < 1 and small, there is just E 0 in the curve C, so the index of the closed curve C is as follows: where I 0 is the index of E 0 . E 0 is a stable node, so I 0 = 1 and: Now we assume that R 0 is close enough to 1 and that two other fixed points X � exist in the C. Therefore, I C is expressed as follows: where I 1 and I 2 are the indices of the other fixed points X � . Thus I C can be considered as a function of R 0 . Now we claim that the I C (R 0 ) is a continuous function. I C (R 0 ) can be rewritten as an integral in the complex plane as follows: whereC is f(C)(C is considered as a closed curve in C and f as a function C ! CÞ. As a result, we can obtain:  We know that I 1,2 = ±1 and I C (R 0 ) is continuous, so: As a result, one of the fixed points X � must be a saddle point (I = −1) and the other one must be a stable or an unstable node (I = 1), but d > 1 4 and λ 1 + λ 2 < 0, so it must be a stable node. Now when R 0 > 1, E 0 is an unstable fixed point (saddle point) and there is always an other fixed point X � 1 in the first quadrant and a fixed point X � 2 near E 0 in the semicircle(X � 2 is not a permissible fixed point because it exists in the fourth quadrant). Therefore, I 0 = −1 and I C (R 0 ) is continuous: As a result, X � 1;2 are nodes. However, it is similar to the case R 0 < 1 and the permissible node must be stable, so X � 1 is a stable fixed point. Therefore, there is just a stable fixed point (E 0 ) for R 0 < R ð1Þ 0 and two stable fixed points (E 0 ; X � 1 ) and an unstable fixed point (X � 2 ) for R ð1Þ 0 < R 0 < 1 and a stable fixed point X � 1 and an unstable fixed point E 0 for R 0 > 1. There must hence be a backward bifurcation when R 0 ¼ R ð1Þ 0 < 1: As an example, this theorem is illustrated in the plots of Fig 4.  remark 0.3. The saddle node bifurcation in backward bifurcation occurs when the quadratic Eq (10) has just one root, a.e: Now when β is changing and other parameters are constant, β can be obtained form the above equation. For example, β = 0.892857 in the simulation of previous theorem, Fig 4. Hopf bifurcation remark 0.4. In the following theorem, R 0 is considered as a function of β, and the other parameters are constant. In other words, as R 0 / β, there is no difference between changing of R 0 and β. It is advisable to define β 0 as follows: Theorem 0.8 (Hopf bifurcation). Suppose bð2dþ1Þ 2 < A and d ðdþ1Þb < b. Define Δ and β max as follows: By considering the above assumptions, the system (6) undergoes a Hopf bifurcation for some β 2 (β 0 , β max ) or equally for some R 0 2 1; A Proof. First, we consider λ 1 + λ 2 (I). We have shown: At First, we fix β 2 [β 0 , β max ]. Now when I = 0, λ 1 + λ 2 = −d < 0 and when I = b:  Tables 1 and 2; R 0 is defined in the Eq (9). https://doi.org/10.1371/journal.pone.0276969.g004 we will have: By considering the above assumptions, for each β 2 [β 0 , β max ], we have λ 1 + λ 2 (0)<0 and λ 1 + λ 2 (b)>0. Therefore, and according to the intermediate value theorem, for each β 2 [β 0 , β max ], there is some I 1 2 (0, b) such that λ 1 + λ 2 (I 1 ) = 0. However, if we consider I 1 as the intersection point of the function Ib ðIþbÞ 2 and the line βI + d, and by considering the behavior of the function and line, there is just one I 1 for each β 2 [β 0 , β max ] and I 1 (β) is continuous and 0 < I 1 (β)<b and λ 1 + λ 2 (I 1 (β)) = 0 when β 2 [β 0 , β max ]. Now we consider the quadratic Eq (10). When d ðdþ1Þb < b and R 0 = 1 or equally β = β 0 , there is just one permissible solution: I = 0; and when R 0 > 1, there is always just one permissible root that is a continuous function of R 0 or other parameters like β (I(β)). Now we consider the root of the quadratic Eq (10) for β = β max . In this case, the quadratic Eq (10) is expressed as follows: We then consider P M (b) and the condition bð2dþ1Þ 2 < A: Obviously, P M (I) has a root I in the interval (b, +1). In other words, by assuming the conditions d As a result, there is a β 1 2 (β 0 , β max ) such that f(β 1 ) = 0 or equally I 1 (β 1 ) = I � (β 1 ) and By considering the differentiability of I(β) and λ 1 + λ 2 (I) and the intermediate value theorem, it is possible to find a β 1 with the following properties: We now consider λ 1 λ 2 (I). We have proved that λ 1 λ 2 (I) is positive when I 2 (0, +1) and d ðdþ1Þb < b, so: Obviously, λ 1,2 are continuous functions of I and β, and I is a continuous function of β, so λ 1,2 (β) is continuous with values in the complex plane C. Consequently, we can select a �, as defined in the above arguments, with the following properties: As l 2 ¼ � l 1 , the above results can be summarized and rewritten as follows: As a result, there must be a Hopf bifurcation when β = β 1 2 (β 0 , β max ) or equally As an example, this theorem can be observe in the plot of Fig 5. It is straightforward to examine that these parameters satisfy the conditions of this theorem, and β 0 = .04 and β max �.06154. Obviously, there are two Hopf bifurcations in this case. These bifurcations occur as β 1 � 0.040968 and β 2 � 0.078668 or equally as R ð1Þ 0 � 1:0242 and R ð2Þ 0 � 1:9667. The case β 1 was predicted in our theorem.

Numerical analysis
In this section, we concentrate on a range of the parameters where there is a Hopf bifurcation, and we subsequently analyze the time evolution of the model. For instance, we can consider the case in Fig 5, where β changes in the interval (.04, .09) with obviously one limit cycle in this situation for β 2 (0.040968, 0.078668). Now we choose two different values for β, before (Fig 6) and after (Fig 7) β 1 � 0.040968, with two different initial conditions for each β, and we sketch the stream plot S − I and S, I, H, R and the fluxes as functions of t, where the fluxes are defined as follows: In Fig 6, β is less than 0.040968, so there is no limit cycle in this condition and the curves in the stream plots (panels (A, A 0 )) approach a fixed point; therefore, the plots of S, I, H, R and the fluxes have a limit when t ! 1. However, in Fig 7, β > 0.040968 and there is a limit cycle and the curves in the stream plots (panels (A, A 0 )) approach this limit cycle. As a result, the plots of S, I, H, R and the fluxes behave periodically when t is large enough.

Discussion & conclusion
In our study, we modified the standard SIR model and divided the infected individuals who survived in the block R and in a newly proposed block H in order to model the number of individuals hospitalized in an ICU with their hospitalization rates defined respectively as nI and bI Iþb . Our model can be considered as a specific version of the model developed by Zhang et al.  Tables 1 and 2; R 0 is defined in the Eq (9).
The most crucial finding of our study is that we demonstrated the existence of Hopf bifurcation and the limit cycle. The existence of the limit cycle, in particular, illustrates that a disease can survive in a community and that this is a potentially critical situation for governments in controlling the spread of a disease.
As an application of our model, we imagine the following case based on the existence of Hopf bifurcation. Suppose that we are in a range of the parameters where a supercritical Hopf bifurcation can occur when β = β 0 , resembling the case considered in the numerical analysis section. Once a given disease spreads, governments usually utilize interventions in order to reduce contact among people, for instance, by implementing quarantine measures. As a result, β is small and we can suppose that β < β 0 where there is no limit cycle. However, after a while, it is possible that β increases and β > β 0 . This event is highly probable to occur, for example, when people ignore interventions of governments and quarantine requirements. There exists, in this situation, a limit cycle and the disease can survive.
A limitation of our study is that we did not analyze the period of limit cycle and the time evolution of the fluxes, which were defined in the Eq (15), which can be researched by numerical analysis more precisely and more convenient instead of real analysis. It is worth mentioning that the fluxes are crucial since they are practically observable in hospitals and can be deployed as a sort of warning signal. Moreover, these fluxes can help estimate the other parameters of the disease spread dynamics. Finally, the proved theorems and the introduced fluxes provide suitable materials for future numerical researches such as global stability, sensitivity, optimal control, which can involve real situations and data of real infection.